#### 这个是学习b站视频 目前李龙源导师还没给ppt  2024年7月12日15:04:11 ####
 ## https://www.bilibili.com/video/BV1nc411g7G3  根据他的可成进行下载测试学习 ## 
library(tidyverse)
# R语言读取gtf文件
library("rtracklayer")
chooseBioCmirror()
# 默认选择国内
4
chooseCRANmirror()
15
# a= c("a","a","ds","a","ds","a","ds","a","ds")
# 
# duplicated(a)
# 
# c<- unique(a)

try({
  setwd('TCGA/')
})

if (!requireNamespace("BiocManager", quietly = TRUE)){
  install.packages("BiocManager")
}

if (!requireNamespace("reshape2", quietly = TRUE)){
  install.packages("reshape2")
}
library(reshape2)

if (!requireNamespace("rtracklayer", quietly = TRUE)){
  BiocManager::install("GSVA")
 # BiocManager::install("rtracklayer")
}

if (!requireNamespace("estimate", quietly = TRUE)){
  BiocManager::install("estimate")
}

# rfore = 'http://r-forge.r-project.org'
# install.packages("estimate",repos = rfore,dependencies = TRUE)

library(BiocManager)
library(estimate)
# chooseBioCmirror()
# BiocManager::install()


# 自动下载TCGAbiolinks 会根据依赖全部内容 不用一个个进行下载 方便
# BiocManager::install("TCGAbiolinks")

library(TCGAbiolinks)
# 获取全部数据
# projects <- TCGAbiolinks::getGDCprojects()$project_id ##获取癌症名字
# # 过滤掉不是以TCGA 开头的数据库
# projects <- projects[grepl('^TCGA', projects, perl=TRUE)]

# 准备下载什么 后面参数暂时使用这个为
# variable = 'TCGA-LUAD'
# str2 <- "_download.rda"
# downFile = paste0(variable, str2)
# downFile
# query <- GDCquery(
#   project = variable,
#   data.category = "Transcriptome Profiling",
#   access = "open",
#   data.type = "Gene Expression Quantification"
# )
# # query <- GDCquery(
# #   project = variable,
# #   data.category = "Simple Nucleotide Variation",
# #   access = "open",
# #   data.type = "Masked Somatic Mutation"
# # )
# 
# #存放的目录  可能会下载很久 
# GDCdownload(query,directory = "GDCdata")
# #存储文件到里面保存成我们需要的数据格式
# expquery2 <- GDCprepare(query = query)
# save(expquery2,file = 'luad.gdc_2022.rda')

load('luad.gdc_2022.rda')

# 
# TCGALUAD load("F:/rproject/r-language/lilongyuan2/TCGA/TCGA-LUAD_download.rda")

gene_annotation_2022 <- ""
if(file.exists("gene_annotation_2022.csv")){
  print("本地存在直接获取")
  gene_annotation_2022 <- read.csv('gene_annotation_2022.CSV')
  gene_annotation_2022<- gene_annotation_2022[,-1]
}else{
  print("不存在进行解析生成到本地文件")
  gene_annotationfROM <- import('gencode.v22.annotation.gtf.gz') 
  gene_annotation <- as.data.frame(gene_annotationfROM)#将文件转换为数据框格式
  gene_annotation <- gene_annotation [gene_annotation$type == 'gene', ] # 筛选为gene的类型
  gene_annotation <- dplyr::select(gene_annotation, gene_id, gene_name,  gene_type)
  colnames(gene_annotation) <- c("ENSEMBL", "symbol","type")
  gene_annotation_2022 <- as.data.frame(gene_annotation)
  write.csv(gene_annotation_2022,'gene_annotation_2022.CSV')
}
# tcga数据库表达普有好几种格式 counts(只用于做差异性分析使用) 形式  tpms 形式(基因的表达普)
# 除了差异性分析都用tpms 做
# fpkm 格式 
# table(gene_annotation_2022$type)
#每一组基因提取都是一样的 方法和字段 
counts1 <- expquery2@assays@data@listData[["unstranded"]]
colnames(counts1) <- expquery2@colData@rownames
rownames(counts1) <- expquery2@rowRanges@ranges@NAMES
counts1<-as.data.frame(counts1)

counts2 <- counts1
# 添加列
counts1 <- rownames_to_column(counts1,var = 'ENSEMBL')
# 链接
counts1 <- inner_join(counts1,gene_annotation_2022,by = 'ENSEMBL')
#去重
counts1 <- counts1[!duplicated(counts1$symbol),]

# 先把一行变成一列
library(utils)
library(estimate)
library(tidyverse)
# 今天就先这样  2024年7月13日16:17:27  部分确实文件没有 先看下流程逻辑
setwd('ESTIMATE/Survival_data/')
# filterCommonGenes()
# 生存分析 
# 生存分析 一群人分为两组 查看哪个活得长  需要知道肿瘤患者的存活时间
# 生存分析 
# os 下载地址 
#https://xenabrowser.net/datapages/?dataset=survival%2FLUAD_survival.txt&host=https%3A%2F%2Ftcga.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443

survivalData <- os
survivalData<- survivalData[,c(1,3,4)]
survivalData$name <- paste0(survivalData$sample,"A")
table(substring(survivalData$name,14,16))
# name 给行名 第一次使用
rownames(survivalData) <- survivalData$name
survivalData<- survivalData[,c(2,3)]
# 看到 https://www.bilibili.com/video/BV1nc411g7G3 P5 45分钟

# 根据免疫得分进行分组 ,视频讲的是根据一列的值取个中位数 
#高于的就是免疫得分高的样本组 低的就是低样本组
# 看了三节课 最核心的是整理数据进行什么类似的生存曲线什么的只需要调参就可以了

# 下一阶段课程 https://www.bilibili.com/video/BV1nc411g7G3?p=6 开头

# 最核心的是整理数据
#### day4 ####
setwd('TCGA/clinical/')
load(file = '../luad.gdc_2022.rda')

clinical <- as.data.frame(expquery2@colData) %>% .[!duplicated(.$sample),]

# 箱线图

# https://www.bilibili.com/video/BV1nc411g7G3?p=7 2024年7月19日09:06:08 新学习

data(mtcars)
mtcarsda<- mtcars
mtcarsda<- mtcarsda[c(1:3),c(1:3)]

#测试melt
colnames(mtcarsda)
mtcarsdacarb <- melt(mtcarsda,id.vars = 'disp')
# 




# https://www.bilibili.com/video/BV1nc411g7G3?p=8 2024年7月19日13:57:02 学习到18分钟
